Heterogeneous slab thermal dehydration driving warm subduction zone earthquakes

Changing thermal regime is one of the key mechanisms driving seismogenic behaviors at cold megathrusts, but it is difficult to interpret warm subduction zones such as Vanuatu for the temperatures are higher than that accommodates shallow brittle failures. We construct a 3-D thermomechanical model to clarify the thermal structure that controls tectonic seismicity in Vanuatu and predict a warm circumstance associated with abundant seismicity. Results reveal a heterogeneous slab ranging from 300 °C to over 900 °C from the Moho to subvolcanic depth. The subduction seismicity corresponds well to the plate interface where dynamic thermal dehydration is focused. The transformation from hydrated basalts to eclogites along the slab facilitates the occurrence of intense earthquakes and slips. Multistage mineralogical metamorphism affects the dynamic stability of megathrusts and favors the generation of active interplate large events. Therefore, slab thermal dehydration plays a greater role than slab temperature condition in influencing the subduction earthquake distribution in warm subduction systems.


Methods and model settings
Based on the thermomechanical model Stag3d (Tackley and Xie, 2003) and the finite difference method (FDM), we investigate the in situ along-strike slab thermal variation in the Vanuatu arc megathrust area and the petrological metamorphic changes in the subducting oceanic lithosphere.An anelastic liquid approximation and the equations of conservation of mass, momentum, and energy are used in this study: (   +  ⋅ ) =  2  + () 2 +   +   , where  is the pressure deviation from the hydrostatic pressure,  0 is the reference thermal expansivity, ij  ( i , j =1, 2, 3) is the stress tensor, and   is the Kronecker delta.
The energy equation includes the advection term, thermal diffusion term, viscous dissipation term, adiabatic heating term, and radioactive heating term.The main model parameters are listed in Tables S3 and S4.
The viscous flow law for wet olivine (Burkett and Billen, 2010) following laboratory experiments is included in this study.The deformation of olivine occurs by both diffusion creep (  ) and dislocation creep (  ), where each mechanism accommodates a portion of the total strain rate (Hirth and Kohlstedt, 2003): The composite upper mantle viscosity under a given condition is as follows: where ηdf and ηds are the diffusion creep and dislocation creep viscosities for olivine, respectively.The general viscosity law is as follows: Here, Ė = ( is the square root of the second invariant of the strain rate tensor (e.g., Ranalli, 1995), Ta is the absolute temperature, R is the gas constant, and Pl is the lithostatic pressure.The values of the model parameters for the diffusion and dislocation creep of olivine are tabulated in Table S1.The trenchward thermal structure follows the Global Depth and Heat (GDH1) model (Stein and Stein, 1992; Grose and Afonso, 2013): The time-dependent thermal boundary condition includes the initialization time.
(,   ) is the temperature at depth z and plate age toc along the trench, Tm is the lithospheric basal temperature,  is the thermal diffusivity, and d0 is the minimum depth of the adiabatic heating layer.(Yoshii, 1975).The incoming oceanic plate comprises a MORB layer at the top with a thickness of 7 km underlain by an ultramafic rock layer (Hacker et al., 2003).The model dimensions are 1600 km×800 km×400 km (along-arc length×across-arc length×depth) and 80×80×100 grids.The temperature boundary condition agrees with the plate cooling model (Grose and Afonoso, 2013).
The bottom of the slab and the vertical planes are prescribed as adiabatic and permeable, and the top surface is set to be a fixed temperature (0°C) and rigid.The subduction velocities inside a prescribed 3-D constrained volume of the oceanic lithosphere are given based on the kinematic plate subduction modeling method (Ji et (, ) = We specified the subduction time to be at least 20 Myr to ensure that the model reaches a steady thermal state with a temperature variation <10℃ over time with a lapse time of ≥ 5 Myr.We tested the resolution and found that the temperature variance was <1% with a maximum temperature variance <1.9% between 80 × 80 × 100 and 96 × 96 × 100 meshes.We performed sensitivity tests to investigate the robustness of our modeling results and varied the mantle viscosity from 1.0×10 19 Pa s to 1.0×10 21  Using the temperature and pressure provided by the numerical simulation, we estimated each facies domain and the corresponding water content (wt%) at every grid.Through the interpolation method, we obtained the intraslab water content distribution (wt%) at various depths.
To calculate the intraslab slab dehydration rate, the result of the 3-D water content at every grid was utilized.Slab dehydration (wt%/km) indicates the change in rock saturation water content (wt%) via a distance (km) in the subduction direction between neighboring grids.We prescribed the slab to be divided into >70 layers according to the mesh number, with the layer surface parallel to the plate interface.Then, the water content at a point on a layer surface was calculated from the water content at the grids surrounding this point using an interpolation method.Thus, the water content at each point on every layer surface was obtained.Next, the difference in the water content value between two neighboring points in the subduction direction was divided by the point distance (which can be derived using the horizontal distance and vertical distance), and then the slab dehydration at each point on a layer could be obtained.Based on these steps, the intraslab dehydration distribution was calculated well.Slab minerals vary in saturation at each subduction stage.Thus, slab dehydration (wt%/km) reflects the efficiency of fluid production by crystalline breakdown and fluid relaxation during subduction along the slab geometry.
The surface heat flow and Curie depth are correlated, following a theoretical thermal conduction relationship (e.g., Li et al., 2011Li et al., , 2013)): where Qs is the surface heat flow, Tc is the Curie temperature at the Curie depth Zb, T0 is the temperature at the surface elevation Zs, K is the average thermal conductivity of the magnetic layer (Table S4), H0 is the heat production rate at the surface, and hr is the Radioactive heat generation rate in mantle 2.245×10 -13a W•m -3   Standard specific heat at constant pressure 1200 a J•kg -

1 2[
( + ,  + ) − ( + , ) + ( − , ) − ( − , is the subduction velocity, and Δ is the interval between two neighboring nodes along the axes.  ,   , and   indicate the model lengths along the x, y, and z axes, respectively.Surface heat flow observations from the Global Heat Flow Database (Pollack et al., 1993) and heat flow values from Curie point depth estimates (Li et al., 2017) are used to constrain the model (Fig. S1).Our model follows the trench temperature boundary of the plate cooling model (Stein and Stein, 1992; Grose and Afonso, 2013).
Pa s and the mantle density from 3250 kg/m 3 to 3350 kg/m 3 .We present the benchmark model results as deviations from the reference models (ΔT and ΔH2O) and show these results at different depth levels within the oceanic slab.The tests show that mantle density variations (±50 kg/m3) induce small temperature variations of <10°C at depth.Ultramafic mantle rocks such as harzburgite (olivine + orthopyroxene) represent the dominant rock type in mantle wedges and the uppermost oceanic mantle, and depleted lherzolite (olivine + orthopyroxene + clinopyroxene) is considered subordinate(Hacker et al., 2003).Seismological studies support the hypothesis that harzburgite represents the principal rock type in the upper mantle.The observed P-wave speeds fromWhite et al. (1992) for the oceanic lower crust and mantle compared with P-wave speeds for various rocks at 200 MPa(Hacker et al., 2003) indicate that most oceanic uppermost mantle (suboceanic mantle) velocity measurements are best explained in terms of spinel harzburgite mantle composition.Due to the reasons above, in our petrological modeling approach, harzburgite is assumed to be the dominant ultramafic rock.We established a P-T-wt%-facies database according toOmori et al. (2009) (MORB) andHacker et al. (2003) with a P-T grid interval of 0.04 GPa (1.2 km) and 5°C.The temperature and pressure at every P-T grid point were calculated from our 3D thermal model.The pressure (GPa) at every grid point was obtained by converting its depths (km) through the preliminary reference Earth model (PREM) parameters.
characteristic drop-off of heat production.The equation shows a nonlinear inverse relationship between heat flow and Zb.For the oceanic lithosphere, we assume H0 = 1.37 μW/m 3 , hr = 5.0 km, Tc = 550 °C, T0 = 5 °C, and Zs = 4 km.For continents, we take H0 = 2.0 μW/m 3 , hr = 10.0 km, and Zs = -1 km to account for the larger radioactive contribution(Li et al., 2017).High heat flow measurements tend to be correlated with small Curie depths, and vice versa.The continental data can be best fitted with the theoretical curve of an average thermal conductivity K of ~2.5 W/m°C, and most oceanic data are best fitted with an average K = ~2.0W/m°C(Li et al., 2013).These conductivities are compatible with those of granite and basalt (e.g.,Turcotte and Schubert, 2002).Synthetic modeling suggests that the largest error in the estimated Curie depths using the linearized centroid method does not reach 35%, and the uncertainty of the surface heat flow is expected to be <20 mW/m 2 due to the selected fractal exponent and wavenumber bands for linear regressions and the observed surface heat flow in plate convergence zones (e.g.,Li et al., 2013Li et al., , 2017)).The model settings combine effective simulation methods used in active convergence zones, such as reconstructions of oceanic-continental cold subduction in eastern Japan(Ji et al., 2017b(Ji et al., , 2017c)), Hikurangi(Suenaga et al., 2018), and Sumatra(Ji et al., 2021) and warm subduction in southwestern Japan(Ji et al., 2016; Ji and   Yoshioka, 2017), Ryukyu(Suenaga et al., 2021), and Cascadia(Ji et al., 2017a).These studies on 3-D thermomechanical models are focused on subducting plates with topographically changing dip angles.

Fig
Fig. S1 (a) Spatial distribution of the surface heat flow calculated in this study.The observed and calculated surface heat flows are compared along trench-normal profiles "b" to "i" in panels (b) to (i).(b) Observed and calculated heat flows along profile b (y = 180 km) in (a).Purple diamonds

Fig. S2
Fig. S2 Calculated thermal state of the incoming plate at various intraslab depths (measured vertically downward from the slab interface).The red cones indicate active arc volcanoes (Siebert et al., 2010).The colored spheres indicate the intraslab seismic events from January 2000 through December 2010 (IRIS; Trabant et al., 2012), with the color indicating the earthquake magnitude, corresponding to the "Magnitude" color scale at the bottom of the figure, and the size mimicking the rupture dimension.(a) The interface.(b) 8 km below the interface.(c) 16 km below the interface.(d) 24 km below the interface.
Fig. S3(a) Profile a